stanEx3
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}non linear regression with one explanatory variable
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex8-0.stan
single linear regression
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -275175
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 36 -38.1623 0.00367238 0.000290015 1 1 86
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -38.16
## b0 7.96
## b1 1.80
## s 4.09
## m0[1] 33.49
## m0[2] 37.42
## m0[3] 25.27
## m0[4] 43.10
## m0[5] 36.89
## m0[6] 25.44
## m0[7] 41.35
## m0[8] 16.79
## m0[9] 29.90
## m0[10] 19.98
## m0[11] 25.63
## m0[12] 29.00
## m0[13] 38.77
## m0[14] 37.31
## m0[15] 39.65
## m0[16] 32.01
## m0[17] 34.18
## m0[18] 13.78
## m0[19] 22.34
## m0[20] 22.12
## y0[1] 30.23
## y0[2] 37.79
## y0[3] 32.51
## y0[4] 41.73
## y0[5] 44.44
## y0[6] 20.67
## y0[7] 38.04
## y0[8] 16.54
## y0[9] 31.74
## y0[10] 19.97
## y0[11] 24.49
## y0[12] 35.25
## y0[13] 28.08
## y0[14] 27.72
## y0[15] 44.25
## y0[16] 31.67
## y0[17] 38.49
## y0[18] 14.84
## y0[19] 10.39
## y0[20] 22.31
## m1[1] 13.78
## m1[2] 17.04
## m1[3] 20.29
## m1[4] 23.55
## m1[5] 26.81
## m1[6] 30.07
## m1[7] 33.33
## m1[8] 36.58
## m1[9] 39.84
## m1[10] 43.10
## y1[1] 9.63
## y1[2] 16.29
## y1[3] 22.19
## y1[4] 22.72
## y1[5] 27.07
## y1[6] 30.94
## y1[7] 26.22
## y1[8] 31.60
## y1[9] 42.91
## y1[10] 42.25
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.26 -37.96 1.22 1.04 -40.65 -36.95 1.01 502 557
## b0 8.36 8.22 2.87 2.83 3.74 13.12 1.01 317 316
## b1 1.78 1.78 0.22 0.21 1.40 2.13 1.01 348 396
## s 4.59 4.50 0.81 0.74 3.50 6.11 1.00 945 911
## m0[1] 33.50 33.51 1.10 1.04 31.71 35.30 1.01 2024 1283
## m0[2] 37.37 37.37 1.34 1.28 35.17 39.51 1.01 1403 1221
## m0[3] 25.41 25.43 1.18 1.14 23.51 27.33 1.01 510 680
## m0[4] 42.97 42.96 1.86 1.75 39.86 45.97 1.00 713 975
## m0[5] 36.85 36.87 1.30 1.24 34.72 38.94 1.01 1509 1205
## m0[6] 25.58 25.59 1.17 1.14 23.69 27.48 1.00 522 695
## m0[7] 41.25 41.23 1.69 1.61 38.46 43.93 1.00 818 976
## m0[8] 17.05 16.99 1.91 1.89 13.94 20.24 1.01 333 357
## m0[9] 29.97 29.97 1.03 0.99 28.29 31.68 1.01 1300 1196
## m0[10] 20.20 20.17 1.60 1.58 17.62 22.83 1.01 357 398
## m0[11] 25.76 25.78 1.16 1.13 23.90 27.64 1.00 536 704
## m0[12] 29.08 29.08 1.04 1.00 27.40 30.78 1.01 1055 1156
## m0[13] 38.71 38.71 1.46 1.38 36.31 41.02 1.00 1180 1120
## m0[14] 37.26 37.26 1.34 1.27 35.08 39.39 1.01 1424 1172
## m0[15] 39.57 39.57 1.53 1.46 37.04 42.01 1.00 1009 1121
## m0[16] 32.04 32.04 1.05 0.99 30.35 33.75 1.01 1872 1267
## m0[17] 34.19 34.18 1.13 1.06 32.33 36.05 1.01 1988 1286
## m0[18] 14.09 13.98 2.23 2.19 10.46 17.74 1.01 323 364
## m0[19] 22.52 22.50 1.40 1.38 20.28 24.79 1.01 396 445
## m0[20] 22.30 22.27 1.41 1.40 20.03 24.62 1.01 391 445
## y0[1] 33.65 33.62 4.83 4.58 25.71 41.79 1.00 1877 1727
## y0[2] 37.27 37.32 4.82 4.64 29.31 45.20 1.00 2025 1864
## y0[3] 25.39 25.46 4.76 4.61 17.28 32.85 1.00 1616 1816
## y0[4] 43.04 43.15 5.00 4.71 34.87 51.06 1.00 1836 1970
## y0[5] 37.06 37.09 4.77 4.74 29.41 44.75 1.00 1965 1835
## y0[6] 25.75 25.81 4.77 4.71 17.94 33.54 1.00 1833 1947
## y0[7] 41.21 41.16 5.05 4.93 32.96 49.55 1.00 1689 1982
## y0[8] 16.94 17.02 5.00 4.86 8.80 24.96 1.00 1343 1548
## y0[9] 30.06 30.16 4.67 4.67 22.39 37.84 1.00 2204 1976
## y0[10] 20.23 20.35 4.91 4.60 11.60 28.22 1.00 1341 1648
## y0[11] 25.66 25.65 4.79 4.82 17.88 33.49 1.00 1906 1866
## y0[12] 29.02 29.09 4.60 4.65 21.33 36.48 1.00 1963 1743
## y0[13] 38.79 38.54 4.88 4.50 31.21 46.90 1.00 1788 2001
## y0[14] 37.17 37.27 4.72 4.54 29.27 44.77 1.00 1983 2124
## y0[15] 39.51 39.61 4.87 4.73 31.25 47.33 1.00 1794 1876
## y0[16] 32.17 32.18 4.84 4.73 24.18 39.96 1.00 1951 1921
## y0[17] 34.23 34.17 4.90 4.77 26.04 42.18 1.00 2098 1878
## y0[18] 14.03 14.01 5.15 4.79 5.57 22.73 1.00 1001 1218
## y0[19] 22.68 22.81 4.86 4.57 14.75 30.60 1.00 1377 1521
## y0[20] 22.31 22.27 4.77 4.60 14.82 30.09 1.00 1620 1743
## m1[1] 14.09 13.98 2.23 2.19 10.46 17.74 1.01 323 364
## m1[2] 17.29 17.24 1.89 1.87 14.23 20.43 1.01 334 371
## m1[3] 20.50 20.47 1.57 1.55 17.98 23.10 1.01 361 417
## m1[4] 23.71 23.70 1.30 1.25 21.62 25.83 1.01 430 485
## m1[5] 26.92 26.92 1.10 1.10 25.14 28.71 1.00 644 879
## m1[6] 30.13 30.13 1.03 0.99 28.46 31.84 1.01 1352 1174
## m1[7] 33.34 33.34 1.09 1.03 31.56 35.12 1.01 2020 1265
## m1[8] 36.55 36.56 1.28 1.22 34.42 38.61 1.01 1579 1259
## m1[9] 39.76 39.75 1.55 1.48 37.23 42.22 1.00 979 1121
## m1[10] 42.97 42.96 1.86 1.75 39.86 45.97 1.00 713 975
## y1[1] 14.18 14.13 5.13 4.95 5.80 22.44 1.00 1036 1540
## y1[2] 17.32 17.44 5.12 4.83 8.75 25.43 1.00 1160 1717
## y1[3] 20.57 20.63 4.83 4.67 12.64 28.35 1.00 1703 1752
## y1[4] 23.70 23.81 4.79 4.60 15.65 31.42 1.00 1669 1799
## y1[5] 26.75 26.72 4.84 4.69 18.99 34.83 1.00 1998 1808
## y1[6] 30.12 30.05 4.77 4.70 22.49 38.06 1.00 1811 1931
## y1[7] 33.47 33.53 4.71 4.64 25.76 40.97 1.00 2046 1808
## y1[8] 36.50 36.49 4.95 4.62 28.28 44.60 1.00 2008 1922
## y1[9] 39.80 39.82 4.87 4.52 31.71 47.78 1.00 1886 1972
## y1[10] 42.84 42.74 5.02 4.93 34.84 51.31 1.00 1820 1738
quadratic regression
y=b0+b2(x-b1)**2
ex8-1.stan
quadratic regression
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real b2;
real<lower=0> s;
}
model{
y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b2*(x[i]-b1)^2;
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b2*(x1[i]-b1)^2;
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)## Initial log joint probability = -574452
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpgKfISJ/model-14463151a723.stan', line 15, column 2 to column 29)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpgKfISJ/model-14463151a723.stan', line 15, column 2 to column 29)
## Error evaluating model log probability: Non-finite gradient.
## 72 -8.59084 0.000254104 0.00246273 1 1 152
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -8.59
## b0 5.11
## b1 4.01
## b2 0.52
## s 0.93
## m0[1] 5.64
## m0[2] 10.65
## m0[3] 17.22
## m0[4] 5.14
## m0[5] 7.61
## m0[6] 8.58
## m0[7] 10.02
## m0[8] 14.83
## m0[9] 5.51
## m0[10] 5.12
## m0[11] 5.46
## m0[12] 5.20
## m0[13] 17.65
## m0[14] 9.01
## m0[15] 6.50
## m0[16] 12.46
## m0[17] 5.94
## m0[18] 7.64
## m0[19] 7.53
## m0[20] 7.08
## y0[1] 5.18
## y0[2] 10.08
## y0[3] 18.66
## y0[4] 6.64
## y0[5] 8.46
## y0[6] 8.90
## y0[7] 10.15
## y0[8] 15.35
## y0[9] 5.81
## y0[10] 4.79
## y0[11] 4.38
## y0[12] 3.22
## y0[13] 18.89
## y0[14] 8.11
## y0[15] 6.22
## y0[16] 13.23
## y0[17] 5.67
## y0[18] 5.87
## y0[19] 9.43
## y0[20] 7.78
## m1[1] 8.58
## m1[2] 6.70
## m1[3] 5.55
## m1[4] 5.11
## m1[5] 5.40
## m1[6] 6.40
## m1[7] 8.13
## m1[8] 10.58
## m1[9] 13.75
## m1[10] 17.65
## y1[1] 6.38
## y1[2] 6.04
## y1[3] 5.65
## y1[4] 6.04
## y1[5] 6.12
## y1[6] 6.18
## y1[7] 8.22
## y1[8] 10.12
## y1[9] 13.70
## y1[10] 15.84
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -10.91 -10.52 1.56 1.38 -13.94 -9.07 1.00 616 918
## b0 5.11 5.12 0.37 0.36 4.47 5.70 1.01 459 553
## b1 3.99 4.00 0.15 0.14 3.72 4.22 1.00 1056 870
## b2 0.52 0.51 0.05 0.05 0.43 0.60 1.00 750 1058
## s 1.10 1.06 0.20 0.18 0.82 1.47 1.00 763 950
## m0[1] 5.64 5.64 0.33 0.32 5.10 6.18 1.01 537 748
## m0[2] 10.66 10.65 0.37 0.35 10.07 11.27 1.00 1201 1220
## m0[3] 17.19 17.20 0.63 0.58 16.13 18.20 1.00 1635 1267
## m0[4] 5.16 5.17 0.38 0.37 4.49 5.77 1.01 473 514
## m0[5] 7.58 7.58 0.47 0.45 6.81 8.36 1.00 1186 754
## m0[6] 8.54 8.54 0.58 0.57 7.59 9.48 1.00 1259 694
## m0[7] 10.03 10.03 0.36 0.35 9.44 10.63 1.00 1015 1341
## m0[8] 14.81 14.82 0.50 0.46 13.95 15.64 1.00 2130 1319
## m0[9] 5.53 5.55 0.40 0.39 4.83 6.16 1.01 502 525
## m0[10] 5.14 5.14 0.36 0.35 4.50 5.71 1.01 464 606
## m0[11] 5.48 5.50 0.40 0.38 4.79 6.11 1.01 499 525
## m0[12] 5.23 5.24 0.39 0.37 4.54 5.85 1.01 478 529
## m0[13] 17.61 17.62 0.66 0.61 16.51 18.68 1.00 1543 1233
## m0[14] 9.03 9.03 0.37 0.36 8.42 9.63 1.00 815 918
## m0[15] 6.49 6.49 0.37 0.36 5.87 7.08 1.00 854 686
## m0[16] 12.46 12.45 0.40 0.37 11.78 13.13 1.01 1918 1337
## m0[17] 5.94 5.94 0.34 0.33 5.38 6.48 1.01 619 771
## m0[18] 7.66 7.67 0.38 0.38 7.03 8.27 1.00 649 780
## m0[19] 7.50 7.50 0.46 0.44 6.74 8.27 1.00 1174 720
## m0[20] 7.06 7.07 0.42 0.40 6.36 7.75 1.00 1078 694
## y0[1] 5.62 5.66 1.16 1.09 3.74 7.52 1.00 1635 1730
## y0[2] 10.62 10.66 1.15 1.13 8.64 12.42 1.00 1972 1809
## y0[3] 17.21 17.22 1.33 1.33 15.05 19.33 1.00 1940 1823
## y0[4] 5.20 5.21 1.20 1.12 3.20 7.10 1.00 1597 1449
## y0[5] 7.60 7.65 1.23 1.13 5.52 9.57 1.00 1785 1535
## y0[6] 8.52 8.51 1.24 1.25 6.47 10.57 1.00 1682 1698
## y0[7] 10.04 10.08 1.21 1.15 8.00 12.05 1.00 1851 1772
## y0[8] 14.83 14.81 1.22 1.15 12.79 16.85 1.00 1957 1747
## y0[9] 5.54 5.53 1.17 1.12 3.60 7.41 1.00 1591 1743
## y0[10] 5.14 5.14 1.12 1.05 3.25 6.92 1.00 1473 1775
## y0[11] 5.49 5.52 1.14 1.10 3.61 7.34 1.00 1618 1531
## y0[12] 5.22 5.23 1.16 1.09 3.33 7.12 1.00 1656 1805
## y0[13] 17.62 17.60 1.30 1.22 15.47 19.72 1.00 2038 1909
## y0[14] 9.05 9.06 1.20 1.14 7.10 11.00 1.00 1788 1788
## y0[15] 6.50 6.52 1.15 1.10 4.58 8.37 1.00 2026 1845
## y0[16] 12.45 12.47 1.19 1.12 10.54 14.46 1.00 1993 1601
## y0[17] 5.94 5.92 1.13 1.11 4.06 7.77 1.00 1822 1909
## y0[18] 7.68 7.70 1.17 1.10 5.78 9.59 1.00 1801 1721
## y0[19] 7.54 7.51 1.18 1.16 5.59 9.45 1.00 1641 1912
## y0[20] 7.09 7.09 1.18 1.14 5.20 9.06 1.00 1956 1821
## m1[1] 8.54 8.54 0.58 0.57 7.59 9.48 1.00 1259 694
## m1[2] 6.69 6.69 0.38 0.37 6.04 7.31 1.00 948 682
## m1[3] 5.55 5.55 0.33 0.32 5.01 6.08 1.01 517 734
## m1[4] 5.13 5.13 0.36 0.35 4.49 5.71 1.01 465 589
## m1[5] 5.42 5.44 0.40 0.38 4.73 6.05 1.01 494 531
## m1[6] 6.43 6.45 0.40 0.39 5.74 7.05 1.01 557 604
## m1[7] 8.15 8.16 0.38 0.37 7.52 8.76 1.00 699 675
## m1[8] 10.59 10.59 0.37 0.35 10.00 11.20 1.00 1177 1218
## m1[9] 13.74 13.75 0.45 0.41 12.96 14.47 1.00 2164 1303
## m1[10] 17.61 17.62 0.66 0.61 16.51 18.68 1.00 1543 1233
## y1[1] 8.54 8.54 1.28 1.27 6.45 10.62 1.00 1706 1761
## y1[2] 6.67 6.70 1.19 1.11 4.64 8.54 1.00 1827 1962
## y1[3] 5.56 5.54 1.15 1.09 3.67 7.42 1.00 1598 1520
## y1[4] 5.13 5.16 1.18 1.11 3.25 7.07 1.00 1696 1658
## y1[5] 5.39 5.39 1.16 1.09 3.51 7.31 1.00 1695 1722
## y1[6] 6.38 6.37 1.17 1.14 4.42 8.24 1.00 1746 1580
## y1[7] 8.15 8.16 1.18 1.16 6.16 9.96 1.00 1642 1490
## y1[8] 10.61 10.61 1.17 1.13 8.71 12.53 1.00 1965 1902
## y1[9] 13.71 13.69 1.23 1.18 11.67 15.70 1.00 1986 1863
## y1[10] 17.58 17.58 1.31 1.21 15.51 19.76 1.00 1938 1841
both log regression
log y=b0+b1log x x,y>0 y=exp b0 x**b1
ex8-2.stan
both log regression
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector<lower=0>[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~lognormal(b0+b1*log(x),s);
}
generated quantities{
vector[N] m0;
vector<lower=0>[N] y0;
for(i in 1:N){
m0[i]=b0+b1*log(x[i]);
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector<lower=0>[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*log(x1[i]);
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
qplot(log(x),log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -176.919
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -24.7469 0.000170539 2.69434e-06 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -24.75
## b0 1.00
## b1 1.84
## s 0.83
## m0[1] 2.17
## m0[2] 4.86
## m0[3] 3.00
## m0[4] 2.61
## m0[5] 3.65
## m0[6] 5.07
## m0[7] 4.83
## m0[8] 2.22
## m0[9] -0.85
## m0[10] 1.56
## m0[11] 1.94
## m0[12] 5.21
## m0[13] 0.94
## m0[14] 0.16
## m0[15] 4.94
## m0[16] -0.11
## m0[17] -2.66
## m0[18] 4.51
## m0[19] 4.26
## m0[20] 4.57
## y0[1] 5.03
## y0[2] 114.89
## y0[3] 16.36
## y0[4] 9.55
## y0[5] 42.82
## y0[6] 301.65
## y0[7] 366.15
## y0[8] 9.87
## y0[9] 0.75
## y0[10] 5.68
## y0[11] 2.39
## y0[12] 173.30
## y0[13] 1.52
## y0[14] 0.25
## y0[15] 85.51
## y0[16] 0.53
## y0[17] 0.09
## y0[18] 42.11
## y0[19] 516.22
## y0[20] 89.65
## m1[1] -2.66
## m1[2] 1.36
## m1[3] 2.53
## m1[4] 3.24
## m1[5] 3.75
## m1[6] 4.15
## m1[7] 4.48
## m1[8] 4.76
## m1[9] 5.00
## m1[10] 5.21
## y1[1] 0.08
## y1[2] 2.82
## y1[3] 12.47
## y1[4] 61.17
## y1[5] 70.81
## y1[6] 24.59
## y1[7] 54.41
## y1[8] 335.42
## y1[9] 232.40
## y1[10] 382.55mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.47 -26.12 1.30 1.05 -28.90 -25.10 1.00 904 1011
## b0 1.01 1.01 0.26 0.25 0.59 1.41 1.01 823 1046
## b1 1.83 1.83 0.17 0.17 1.55 2.11 1.00 1133 1191
## s 0.94 0.91 0.17 0.15 0.70 1.27 1.00 1260 1200
## m0[1] 2.18 2.17 0.21 0.21 1.82 2.51 1.00 1030 1081
## m0[2] 4.86 4.85 0.30 0.28 4.35 5.36 1.00 2293 1652
## m0[3] 3.00 3.00 0.21 0.20 2.65 3.35 1.00 1473 1465
## m0[4] 2.62 2.62 0.21 0.19 2.27 2.95 1.00 1243 1297
## m0[5] 3.65 3.65 0.23 0.21 3.25 4.04 1.00 1970 1438
## m0[6] 5.06 5.06 0.31 0.30 4.54 5.59 1.00 2264 1596
## m0[7] 4.83 4.83 0.29 0.28 4.33 5.33 1.00 2298 1652
## m0[8] 2.22 2.22 0.21 0.21 1.87 2.56 1.00 1059 1078
## m0[9] -0.84 -0.84 0.39 0.37 -1.46 -0.21 1.01 821 889
## m0[10] 1.56 1.56 0.23 0.23 1.18 1.92 1.01 878 997
## m0[11] 1.94 1.94 0.22 0.21 1.58 2.29 1.01 957 1110
## m0[12] 5.21 5.20 0.32 0.31 4.66 5.75 1.00 2233 1546
## m0[13] 0.95 0.95 0.26 0.25 0.52 1.36 1.01 820 1002
## m0[14] 0.17 0.17 0.31 0.29 -0.34 0.68 1.01 803 869
## m0[15] 4.93 4.93 0.30 0.29 4.42 5.45 1.00 2285 1681
## m0[16] -0.09 -0.09 0.33 0.31 -0.63 0.45 1.01 805 869
## m0[17] -2.64 -2.64 0.53 0.53 -3.50 -1.75 1.01 872 888
## m0[18] 4.50 4.50 0.27 0.26 4.04 4.97 1.00 2325 1596
## m0[19] 4.26 4.25 0.26 0.24 3.82 4.70 1.00 2296 1651
## m0[20] 4.57 4.57 0.28 0.26 4.10 5.05 1.00 2320 1596
## y0[1] 14.53 8.80 19.98 7.58 1.87 43.99 1.00 1984 1892
## y0[2] 205.52 129.33 254.86 112.09 25.23 630.10 1.00 1802 1800
## y0[3] 32.26 20.22 49.33 16.91 4.18 89.64 1.00 2058 1859
## y0[4] 22.71 13.34 32.73 10.89 2.82 68.87 1.00 1871 1892
## y0[5] 60.24 37.58 78.61 31.25 7.61 173.87 1.00 2165 2100
## y0[6] 273.64 153.51 434.34 135.14 31.94 846.11 1.00 2005 1800
## y0[7] 210.84 123.07 368.38 104.13 25.80 627.61 1.00 1957 1877
## y0[8] 14.64 9.24 20.65 7.90 1.78 44.43 1.00 1948 1778
## y0[9] 0.76 0.42 1.76 0.38 0.08 2.33 1.00 1593 1945
## y0[10] 7.83 4.80 13.55 4.05 0.96 22.34 1.00 1562 1965
## y0[11] 11.47 7.04 18.32 5.94 1.45 32.76 1.00 2032 1769
## y0[12] 327.15 188.81 685.13 163.01 33.92 968.30 1.00 2096 1973
## y0[13] 4.53 2.73 6.93 2.37 0.47 13.75 1.00 1868 1723
## y0[14] 2.06 1.15 9.86 1.01 0.22 5.58 1.00 1884 1828
## y0[15] 229.64 136.08 396.83 116.95 27.88 668.82 1.00 2036 1881
## y0[16] 1.48 0.88 1.97 0.74 0.17 4.86 1.00 1685 1574
## y0[17] 0.15 0.07 0.76 0.07 0.01 0.46 1.00 1649 1821
## y0[18] 160.10 92.88 253.38 82.21 17.97 478.66 1.00 2110 1989
## y0[19] 118.89 69.71 224.96 58.31 13.38 351.97 1.00 2084 1868
## y0[20] 175.03 97.40 310.72 87.05 18.99 555.14 1.00 2126 1611
## m1[1] -2.64 -2.64 0.53 0.53 -3.50 -1.75 1.01 872 888
## m1[2] 1.37 1.37 0.24 0.23 0.98 1.74 1.01 853 998
## m1[3] 2.54 2.54 0.21 0.20 2.19 2.87 1.00 1200 1212
## m1[4] 3.24 3.24 0.22 0.20 2.88 3.60 1.00 1654 1396
## m1[5] 3.75 3.75 0.23 0.22 3.35 4.15 1.00 2051 1487
## m1[6] 4.15 4.15 0.25 0.24 3.72 4.59 1.00 2268 1624
## m1[7] 4.48 4.47 0.27 0.25 4.01 4.94 1.00 2327 1596
## m1[8] 4.75 4.75 0.29 0.27 4.26 5.25 1.00 2309 1651
## m1[9] 4.99 4.99 0.31 0.29 4.47 5.51 1.00 2275 1596
## m1[10] 5.21 5.20 0.32 0.31 4.66 5.75 1.00 2233 1546
## y1[1] 0.13 0.07 0.20 0.06 0.01 0.40 1.00 1434 1436
## y1[2] 6.36 3.95 10.29 3.36 0.75 18.00 1.00 1753 1853
## y1[3] 20.34 12.53 27.11 10.50 2.56 60.89 1.00 2158 1935
## y1[4] 40.41 25.49 53.29 21.06 5.28 121.35 1.00 2118 1894
## y1[5] 72.14 42.34 104.31 36.22 8.46 227.50 1.00 1965 1810
## y1[6] 106.71 63.13 143.03 54.91 12.58 327.53 1.00 2081 1965
## y1[7] 152.74 89.73 356.89 77.32 17.86 448.15 1.00 2249 2055
## y1[8] 184.17 116.68 237.58 98.08 22.92 562.03 1.00 1731 1863
## y1[9] 244.00 145.46 311.55 132.50 25.79 781.64 1.00 2225 1930
## y1[10] 319.77 185.35 611.96 163.76 30.88 922.32 1.00 2012 1930
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
qplot(log(x),log(y),col=I('red'))+
geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=m),xy,col='black')
exponential increasing/decreasing
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)
b1>0 x->Infinity,y->Infinity
b1<0 x->Infinity,y->+0 n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)ex8-3-1.stan
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0*exp(b1*x),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*exp(b1*x[i]);
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*exp(b1*x1[i]);
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)## Initial log joint probability = -1.03578e+08
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 7 -19.3078 0.000234829 0.000118616 1 1 39
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -19.31
## b0 0.00
## b1 -10.00
## s 1.59
## m0[1] 0.00
## m0[2] 0.00
## m0[3] 0.00
## m0[4] 0.00
## m0[5] 0.00
## m0[6] 0.00
## m0[7] 0.00
## m0[8] 0.00
## m0[9] 0.00
## m0[10] 0.00
## m0[11] 0.00
## m0[12] 0.00
## m0[13] 0.00
## m0[14] 0.00
## m0[15] 0.00
## m0[16] 0.00
## m0[17] 0.00
## m0[18] 0.00
## m0[19] 0.00
## m0[20] 0.00
## y0[1] -1.06
## y0[2] 0.29
## y0[3] -0.11
## y0[4] -0.69
## y0[5] -0.80
## y0[6] -0.47
## y0[7] 1.21
## y0[8] 0.46
## y0[9] -1.59
## y0[10] -0.52
## y0[11] -1.15
## y0[12] 1.43
## y0[13] -0.03
## y0[14] 0.92
## y0[15] 1.13
## y0[16] -3.54
## y0[17] -0.07
## y0[18] -0.68
## y0[19] -1.11
## y0[20] -1.20
## m1[1] 0.00
## m1[2] 0.00
## m1[3] 0.00
## m1[4] 0.00
## m1[5] 0.00
## m1[6] 0.00
## m1[7] 0.00
## m1[8] 0.00
## m1[9] 0.00
## m1[10] 0.00
## y1[1] 2.52
## y1[2] -0.90
## y1[3] 1.30
## y1[4] 2.21
## y1[5] 1.79
## y1[6] 1.62
## y1[7] 0.86
## y1[8] -0.35
## y1[9] 2.94
## y1[10] 0.58
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 6.84 7.22 1.40 1.11 4.10 8.35 1.00 456 456
## b0 14.29 12.96 5.44 3.91 8.32 25.30 1.00 356 318
## b1 -2.74 -2.66 0.63 0.57 -3.98 -1.85 1.00 366 321
## s 0.54 0.52 0.10 0.09 0.41 0.71 1.01 491 608
## m0[1] 0.00 0.00 0.01 0.00 0.00 0.01 1.00 371 355
## m0[2] 0.00 0.00 0.01 0.00 0.00 0.01 1.00 371 355
## m0[3] 0.05 0.04 0.05 0.04 0.00 0.15 1.00 377 400
## m0[4] 0.01 0.01 0.02 0.01 0.00 0.05 1.00 374 379
## m0[5] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 371 355
## m0[6] 0.36 0.35 0.17 0.17 0.11 0.68 1.00 392 410
## m0[7] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 370 333
## m0[8] 2.08 2.10 0.27 0.25 1.57 2.49 1.00 689 627
## m0[9] 0.01 0.00 0.01 0.00 0.00 0.03 1.00 373 379
## m0[10] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 371 355
## m0[11] 0.02 0.01 0.02 0.01 0.00 0.06 1.00 374 379
## m0[12] 0.01 0.01 0.02 0.01 0.00 0.05 1.00 374 379
## m0[13] 4.28 4.25 0.43 0.40 3.60 5.03 1.00 603 794
## m0[14] 2.31 2.32 0.26 0.24 1.83 2.71 1.00 846 637
## m0[15] 0.06 0.05 0.06 0.04 0.01 0.18 1.00 377 400
## m0[16] 1.70 1.73 0.28 0.26 1.18 2.13 1.00 542 552
## m0[17] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 371 355
## m0[18] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 370 333
## m0[19] 3.56 3.56 0.30 0.29 3.10 4.06 1.00 1149 1089
## m0[20] 1.32 1.33 0.28 0.27 0.81 1.77 1.00 468 521
## y0[1] -0.02 -0.01 0.55 0.53 -0.92 0.88 1.00 2065 2010
## y0[2] 0.01 0.02 0.54 0.51 -0.88 0.87 1.00 2090 1782
## y0[3] 0.05 0.05 0.55 0.53 -0.86 0.93 1.00 2033 1848
## y0[4] 0.01 -0.02 0.55 0.53 -0.86 0.91 1.00 1959 1857
## y0[5] -0.01 0.00 0.56 0.53 -0.94 0.88 1.00 2117 1972
## y0[6] 0.37 0.35 0.59 0.57 -0.61 1.31 1.00 1583 1738
## y0[7] 0.03 0.03 0.55 0.52 -0.88 0.92 1.00 1763 1703
## y0[8] 2.07 2.08 0.62 0.58 1.01 3.04 1.00 1625 1599
## y0[9] 0.03 0.01 0.55 0.51 -0.86 0.92 1.00 1958 1678
## y0[10] -0.01 -0.01 0.54 0.52 -0.90 0.87 1.00 2007 1883
## y0[11] 0.01 0.02 0.55 0.52 -0.85 0.91 1.00 2014 1892
## y0[12] 0.01 0.03 0.54 0.53 -0.87 0.85 1.00 1804 1694
## y0[13] 4.28 4.27 0.70 0.69 3.18 5.49 1.00 1058 1088
## y0[14] 2.30 2.33 0.61 0.56 1.25 3.30 1.00 1289 1564
## y0[15] 0.06 0.06 0.52 0.51 -0.78 0.90 1.00 1910 1892
## y0[16] 1.70 1.72 0.61 0.56 0.70 2.67 1.00 1296 1528
## y0[17] 0.01 0.01 0.54 0.51 -0.88 0.88 1.00 1986 1877
## y0[18] 0.01 0.01 0.54 0.51 -0.87 0.92 1.00 1976 1973
## y0[19] 3.55 3.56 0.61 0.58 2.53 4.55 1.00 1454 1306
## y0[20] 1.33 1.34 0.62 0.59 0.32 2.30 1.00 1420 1595
## m1[1] 4.28 4.25 0.43 0.40 3.60 5.03 1.00 603 794
## m1[2] 1.15 1.16 0.27 0.27 0.66 1.59 1.00 447 454
## m1[3] 0.33 0.32 0.16 0.16 0.10 0.64 1.00 390 410
## m1[4] 0.10 0.09 0.08 0.06 0.01 0.26 1.00 380 400
## m1[5] 0.03 0.02 0.04 0.02 0.00 0.11 1.00 375 379
## m1[6] 0.01 0.01 0.02 0.01 0.00 0.04 1.00 373 379
## m1[7] 0.00 0.00 0.01 0.00 0.00 0.02 1.00 372 355
## m1[8] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 371 355
## m1[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 370 333
## m1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 370 333
## y1[1] 4.27 4.26 0.68 0.66 3.17 5.40 1.00 1211 1220
## y1[2] 1.18 1.18 0.62 0.59 0.13 2.15 1.00 1180 1120
## y1[3] 0.33 0.32 0.57 0.55 -0.58 1.28 1.00 1563 1505
## y1[4] 0.09 0.10 0.56 0.54 -0.83 1.00 1.00 2113 1447
## y1[5] 0.04 0.03 0.54 0.52 -0.84 0.93 1.00 2012 1794
## y1[6] 0.03 0.04 0.53 0.53 -0.81 0.88 1.00 2085 1892
## y1[7] 0.00 0.00 0.54 0.51 -0.91 0.88 1.00 2073 1706
## y1[8] -0.02 -0.01 0.55 0.55 -0.92 0.87 1.00 2030 1882
## y1[9] 0.00 0.00 0.54 0.51 -0.86 0.87 1.00 2009 1887
## y1[10] -0.01 -0.02 0.54 0.52 -0.89 0.88 1.00 2084 1684
ex8-3-2.stan
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real<lower=-10,upper=10> b1;
real<lower=0> s;
}
model{
y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=log(b0)+b1*x[i];
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=log(b0)+b1*x1[i];
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1167.44
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 36 -14.706 0.000243316 0.000413939 1 1 56
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.71
## b0 10.50
## b1 -2.03
## s 0.50
## m0[1] -6.99
## m0[2] 0.60
## m0[3] -4.02
## m0[4] -6.57
## m0[5] -3.17
## m0[6] 0.44
## m0[7] -4.85
## m0[8] -0.78
## m0[9] 1.80
## m0[10] -1.18
## m0[11] -5.96
## m0[12] 1.56
## m0[13] -6.64
## m0[14] -0.34
## m0[15] -6.65
## m0[16] -3.46
## m0[17] -1.27
## m0[18] 1.37
## m0[19] -5.78
## m0[20] -3.70
## y0[1] 0.00
## y0[2] 2.44
## y0[3] 0.03
## y0[4] 0.00
## y0[5] 0.05
## y0[6] 1.09
## y0[7] 0.01
## y0[8] 0.54
## y0[9] 3.55
## y0[10] 0.43
## y0[11] 0.00
## y0[12] 7.80
## y0[13] 0.00
## y0[14] 0.95
## y0[15] 0.00
## y0[16] 0.04
## y0[17] 0.24
## y0[18] 3.79
## y0[19] 0.00
## y0[20] 0.01
## m1[1] 1.80
## m1[2] 0.82
## m1[3] -0.16
## m1[4] -1.13
## m1[5] -2.11
## m1[6] -3.08
## m1[7] -4.06
## m1[8] -5.04
## m1[9] -6.01
## m1[10] -6.99
## y1[1] 9.00
## y1[2] 1.75
## y1[3] 0.73
## y1[4] 0.21
## y1[5] 0.13
## y1[6] 0.02
## y1[7] 0.04
## y1[8] 0.00
## y1[9] 0.00
## y1[10] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.12 -12.80 1.34 1.16 -15.69 -11.63 1.00 632 894
## b0 11.42 10.94 3.11 2.68 7.38 17.27 1.01 483 529
## b1 -2.04 -2.04 0.09 0.08 -2.19 -1.90 1.00 639 818
## s 0.58 0.57 0.11 0.10 0.43 0.77 1.01 650 840
## m0[1] -7.01 -7.00 0.23 0.22 -7.39 -6.64 1.00 1657 1535
## m0[2] 0.64 0.62 0.20 0.19 0.32 0.99 1.01 502 554
## m0[3] -4.02 -4.02 0.15 0.14 -4.27 -3.78 1.00 2002 1261
## m0[4] -6.58 -6.58 0.21 0.20 -6.94 -6.24 1.00 1776 1614
## m0[5] -3.17 -3.17 0.14 0.13 -3.39 -2.93 1.00 1380 1189
## m0[6] 0.48 0.47 0.19 0.18 0.17 0.83 1.01 505 570
## m0[7] -4.85 -4.85 0.17 0.15 -5.13 -4.59 1.00 2212 1517
## m0[8] -0.76 -0.76 0.16 0.15 -1.00 -0.48 1.00 571 623
## m0[9] 1.84 1.83 0.24 0.23 1.47 2.26 1.01 486 518
## m0[10] -1.16 -1.16 0.15 0.14 -1.40 -0.89 1.00 616 709
## m0[11] -5.97 -5.97 0.20 0.18 -6.30 -5.65 1.00 1968 1616
## m0[12] 1.60 1.59 0.23 0.22 1.24 2.00 1.01 487 518
## m0[13] -6.66 -6.65 0.22 0.20 -7.02 -6.30 1.00 1755 1614
## m0[14] -0.30 -0.31 0.17 0.16 -0.57 0.00 1.00 538 636
## m0[15] -6.67 -6.66 0.22 0.20 -7.02 -6.31 1.00 1753 1614
## m0[16] -3.46 -3.46 0.14 0.13 -3.69 -3.23 1.00 1601 1269
## m0[17] -1.25 -1.26 0.15 0.14 -1.49 -0.98 1.00 629 709
## m0[18] 1.41 1.40 0.22 0.21 1.05 1.81 1.01 488 518
## m0[19] -5.79 -5.79 0.19 0.18 -6.11 -5.49 1.00 2024 1589
## m0[20] -3.69 -3.69 0.14 0.13 -3.93 -3.46 1.00 1777 1284
## y0[1] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2064 1895
## y0[2] 2.29 1.84 1.75 1.02 0.69 5.27 1.00 1316 1799
## y0[3] 0.02 0.02 0.01 0.01 0.01 0.05 1.00 1987 1843
## y0[4] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1957 1792
## y0[5] 0.05 0.04 0.03 0.02 0.02 0.11 1.00 2017 1879
## y0[6] 2.01 1.67 1.49 0.94 0.62 4.52 1.00 1611 1799
## y0[7] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 1931 1919
## y0[8] 0.58 0.48 0.41 0.26 0.18 1.35 1.00 1646 1893
## y0[9] 7.85 6.32 5.81 3.65 2.31 18.10 1.00 1378 1095
## y0[10] 0.38 0.31 0.27 0.18 0.11 0.86 1.00 1789 1680
## y0[11] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1806 1933
## y0[12] 5.97 4.88 4.24 2.69 1.79 13.91 1.00 1445 1516
## y0[13] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1986 1514
## y0[14] 0.91 0.74 0.79 0.43 0.26 2.05 1.00 1748 1448
## y0[15] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2054 1773
## y0[16] 0.04 0.03 0.03 0.02 0.01 0.08 1.00 2130 1838
## y0[17] 0.34 0.28 0.24 0.16 0.10 0.76 1.00 1624 1919
## y0[18] 5.07 4.12 3.86 2.35 1.52 11.78 1.00 1360 1526
## y0[19] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 2051 1816
## y0[20] 0.03 0.03 0.02 0.01 0.01 0.07 1.00 1748 1915
## m1[1] 1.84 1.83 0.24 0.23 1.47 2.26 1.01 486 518
## m1[2] 0.86 0.85 0.21 0.19 0.53 1.22 1.01 497 541
## m1[3] -0.12 -0.13 0.18 0.17 -0.40 0.19 1.00 528 636
## m1[4] -1.11 -1.11 0.16 0.14 -1.35 -0.84 1.00 610 683
## m1[5] -2.09 -2.09 0.14 0.13 -2.31 -1.85 1.00 816 906
## m1[6] -3.07 -3.08 0.14 0.13 -3.30 -2.84 1.00 1313 1179
## m1[7] -4.06 -4.06 0.15 0.14 -4.31 -3.82 1.00 1999 1237
## m1[8] -5.04 -5.04 0.17 0.16 -5.33 -4.77 1.00 2203 1513
## m1[9] -6.02 -6.02 0.20 0.18 -6.35 -5.71 1.00 1953 1616
## m1[10] -7.01 -7.00 0.23 0.22 -7.39 -6.64 1.00 1657 1535
## y1[1] 7.80 6.32 5.96 3.46 2.27 18.72 1.00 1376 1585
## y1[2] 2.80 2.32 1.92 1.33 0.83 6.54 1.00 1575 1681
## y1[3] 1.06 0.87 0.86 0.49 0.31 2.32 1.00 1654 1506
## y1[4] 0.40 0.33 0.26 0.18 0.13 0.88 1.00 1846 1951
## y1[5] 0.15 0.12 0.10 0.07 0.05 0.34 1.00 1801 1891
## y1[6] 0.05 0.05 0.04 0.03 0.02 0.12 1.00 2081 1841
## y1[7] 0.02 0.02 0.02 0.01 0.01 0.05 1.00 2264 1866
## y1[8] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 1725 1929
## y1[9] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1950 1732
## y1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1933 1691
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,log(y),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
growth curve
y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)ex8-4.stan
growth curve
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=10> b1;
real<lower=0> s;
}
model{
y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(1-exp(-b1*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(1-exp(-b1*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)## Initial log joint probability = -3396.36
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -49.0744 0.00344815 0.00157874 1 1 15
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -49.07
## b0 0.00
## b1 1.84
## s 7.05
## m0[1] 0.00
## m0[2] 0.00
## m0[3] 0.00
## m0[4] 0.00
## m0[5] 0.00
## m0[6] 0.00
## m0[7] 0.00
## m0[8] 0.00
## m0[9] 0.00
## m0[10] 0.00
## m0[11] 0.00
## m0[12] 0.00
## m0[13] 0.00
## m0[14] 0.00
## m0[15] 0.00
## m0[16] 0.00
## m0[17] 0.00
## m0[18] 0.00
## m0[19] 0.00
## m0[20] 0.00
## y0[1] -8.04
## y0[2] 10.33
## y0[3] 3.92
## y0[4] 4.07
## y0[5] -4.23
## y0[6] 4.83
## y0[7] 11.15
## y0[8] -1.18
## y0[9] -3.67
## y0[10] -0.31
## y0[11] 5.87
## y0[12] -17.71
## y0[13] -4.30
## y0[14] -0.12
## y0[15] -12.58
## y0[16] 8.36
## y0[17] 1.83
## y0[18] 6.85
## y0[19] 3.86
## y0[20] 2.79
## m1[1] 0.00
## m1[2] 0.00
## m1[3] 0.00
## m1[4] 0.00
## m1[5] 0.00
## m1[6] 0.00
## m1[7] 0.00
## m1[8] 0.00
## m1[9] 0.00
## m1[10] 0.00
## y1[1] -4.59
## y1[2] -10.28
## y1[3] -10.30
## y1[4] 0.40
## y1[5] 0.81
## y1[6] -1.43
## y1[7] -4.01
## y1[8] -3.04
## y1[9] 1.49
## y1[10] 3.87
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.47 -11.10 1.37 1.05 -14.17 -10.03 1.00 642 798
## b0 9.70 9.60 0.86 0.75 8.51 11.26 1.00 737 783
## b1 0.44 0.43 0.10 0.09 0.29 0.61 1.01 621 682
## s 1.20 1.17 0.22 0.19 0.90 1.61 1.00 986 839
## m0[1] 1.60 1.58 0.22 0.21 1.25 1.99 1.00 626 760
## m0[2] 4.78 4.76 0.44 0.42 4.06 5.52 1.00 666 853
## m0[3] 6.54 6.54 0.40 0.39 5.88 7.18 1.00 838 898
## m0[4] 5.05 5.03 0.44 0.42 4.33 5.78 1.00 675 881
## m0[5] 6.69 6.69 0.39 0.38 6.04 7.31 1.00 882 915
## m0[6] 6.74 6.74 0.39 0.38 6.09 7.35 1.00 897 971
## m0[7] 9.16 9.17 0.51 0.49 8.33 10.01 1.00 961 1204
## m0[8] 7.26 7.27 0.35 0.34 6.68 7.83 1.00 1211 1015
## m0[9] 9.05 9.06 0.47 0.45 8.28 9.82 1.00 1052 1286
## m0[10] 9.16 9.17 0.51 0.49 8.33 10.02 1.00 958 1204
## m0[11] 4.72 4.70 0.44 0.42 4.00 5.45 1.00 665 852
## m0[12] 8.09 8.08 0.32 0.31 7.58 8.62 1.00 1990 1195
## m0[13] 8.87 8.88 0.42 0.39 8.19 9.55 1.00 1239 1161
## m0[14] 4.74 4.72 0.44 0.42 4.02 5.48 1.00 665 852
## m0[15] 7.69 7.69 0.33 0.31 7.16 8.23 1.00 1659 1174
## m0[16] 9.19 9.19 0.52 0.51 8.35 10.06 1.00 941 1204
## m0[17] 1.83 1.81 0.25 0.24 1.44 2.27 1.00 627 779
## m0[18] 8.97 8.98 0.45 0.43 8.25 9.71 1.00 1126 1269
## m0[19] 8.25 8.24 0.33 0.31 7.73 8.78 1.00 1999 1184
## m0[20] 1.19 1.17 0.17 0.16 0.92 1.49 1.00 625 759
## y0[1] 1.59 1.57 1.21 1.17 -0.36 3.60 1.00 1956 1645
## y0[2] 4.73 4.74 1.31 1.26 2.57 6.88 1.00 1757 1997
## y0[3] 6.52 6.54 1.30 1.24 4.39 8.70 1.00 2007 2012
## y0[4] 5.05 5.04 1.30 1.26 2.96 7.26 1.00 1862 1773
## y0[5] 6.67 6.64 1.26 1.21 4.68 8.74 1.00 1993 1880
## y0[6] 6.73 6.72 1.30 1.24 4.57 8.92 1.00 1581 1892
## y0[7] 9.15 9.14 1.28 1.26 7.06 11.29 1.00 1751 1777
## y0[8] 7.25 7.25 1.26 1.21 5.15 9.34 1.00 1800 1932
## y0[9] 9.07 9.10 1.31 1.28 6.89 11.19 1.00 1994 1903
## y0[10] 9.16 9.12 1.32 1.24 7.01 11.41 1.00 1647 1808
## y0[11] 4.74 4.74 1.27 1.22 2.68 6.89 1.00 1522 1711
## y0[12] 8.08 8.06 1.29 1.16 6.03 10.23 1.00 1651 1821
## y0[13] 8.84 8.85 1.31 1.24 6.67 10.98 1.00 1881 1657
## y0[14] 4.75 4.75 1.31 1.24 2.60 6.87 1.00 1434 1644
## y0[15] 7.67 7.70 1.30 1.26 5.60 9.74 1.00 1856 1734
## y0[16] 9.17 9.18 1.31 1.28 6.98 11.32 1.00 1865 1822
## y0[17] 1.84 1.81 1.25 1.20 -0.18 3.86 1.00 1971 1820
## y0[18] 8.96 9.00 1.28 1.23 6.77 11.01 1.00 1691 1840
## y0[19] 8.28 8.30 1.28 1.21 6.15 10.38 1.00 2064 2099
## y0[20] 1.20 1.17 1.25 1.20 -0.78 3.18 1.00 2112 1969
## m1[1] 1.19 1.17 0.17 0.16 0.92 1.49 1.00 625 759
## m1[2] 3.60 3.58 0.40 0.38 2.97 4.30 1.00 642 793
## m1[3] 5.31 5.29 0.44 0.43 4.58 6.04 1.00 687 881
## m1[4] 6.52 6.52 0.40 0.39 5.85 7.16 1.00 835 898
## m1[5] 7.38 7.39 0.35 0.33 6.83 7.94 1.00 1330 1134
## m1[6] 8.00 7.99 0.32 0.30 7.48 8.52 1.00 1942 1222
## m1[7] 8.45 8.45 0.35 0.32 7.90 9.00 1.00 1840 1253
## m1[8] 8.78 8.79 0.40 0.38 8.13 9.42 1.00 1354 1287
## m1[9] 9.01 9.02 0.46 0.44 8.27 9.76 1.00 1084 1268
## m1[10] 9.19 9.19 0.52 0.51 8.35 10.06 1.00 941 1204
## y1[1] 1.21 1.21 1.25 1.17 -0.83 3.26 1.00 1711 1840
## y1[2] 3.59 3.60 1.29 1.25 1.42 5.68 1.00 1724 1893
## y1[3] 5.32 5.33 1.31 1.25 3.21 7.45 1.00 1666 1590
## y1[4] 6.52 6.49 1.33 1.25 4.38 8.72 1.00 1618 1747
## y1[5] 7.44 7.44 1.26 1.29 5.36 9.45 1.00 1707 1998
## y1[6] 7.98 7.99 1.28 1.22 5.81 10.02 1.00 1959 1906
## y1[7] 8.48 8.48 1.30 1.29 6.30 10.55 1.00 1646 1353
## y1[8] 8.80 8.81 1.29 1.25 6.70 10.92 1.00 2038 1972
## y1[9] 9.04 9.03 1.28 1.23 6.91 11.11 1.00 1945 1959
## y1[10] 9.19 9.19 1.32 1.27 7.09 11.32 1.00 1607 1314
sigmoid curve
y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)ex8-5.stan
sigmoid curve
data{
int N;
vector[N] x;
int Ym;
array[N] int y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=100> b1;
}
model{
y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
array[N] int y0;
for(i in 1:N){
y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
}
array[N1] int y1;
for(i in 1:N1){
y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
}
}n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-5.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "y0" "y1"seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -59.64 -59.26 1.13 0.75 -62.00 -58.61 1.00 641 796
## b0 3.92 3.93 0.15 0.15 3.66 4.17 1.00 2268 1417
## b1 1.58 1.57 0.24 0.23 1.21 1.99 1.00 564 491
## y0[1] 9.80 10.00 0.46 0.00 9.00 10.00 1.00 1847 NA
## y0[2] 9.98 10.00 0.14 0.00 10.00 10.00 1.00 1781 NA
## y0[3] 9.57 10.00 0.68 0.00 8.00 10.00 1.00 1438 NA
## y0[4] 0.06 0.00 0.26 0.00 0.00 1.00 1.00 1645 1656
## y0[5] 7.86 8.00 1.34 1.48 6.00 10.00 1.00 1179 NA
## y0[6] 0.06 0.00 0.24 0.00 0.00 1.00 1.00 1521 1504
## y0[7] 9.83 10.00 0.43 0.00 9.00 10.00 1.00 1868 NA
## y0[8] 9.80 10.00 0.45 0.00 9.00 10.00 1.00 1765 NA
## y0[9] 0.12 0.00 0.36 0.00 0.00 1.00 1.00 1659 1646
## y0[10] 2.47 2.00 1.47 1.48 0.00 5.00 1.00 1882 1693
## y0[11] 1.15 1.00 1.05 1.48 0.00 3.00 1.00 1790 1978
## y0[12] 2.37 2.00 1.43 1.48 0.00 5.00 1.00 1795 1914
## y0[13] 0.21 0.00 0.47 0.00 0.00 1.00 1.00 1931 1838
## y0[14] 4.97 5.00 1.65 1.48 2.00 8.00 1.00 1958 1776
## y0[15] 8.50 9.00 1.20 1.48 6.00 10.00 1.00 1423 NA
## y0[16] 4.58 5.00 1.62 1.48 2.00 7.00 1.00 2077 1813
## y0[17] 0.69 0.00 0.83 0.00 0.00 2.00 1.00 1862 1705
## y0[18] 5.73 6.00 1.64 1.48 3.00 8.00 1.00 1948 2019
## y0[19] 1.91 2.00 1.30 1.48 0.00 4.00 1.00 1860 1927
## y0[20] 9.37 10.00 0.81 0.00 8.00 10.00 1.00 1863 NA
## y1[1] 0.06 0.00 0.25 0.00 0.00 1.00 1.00 1666 1440
## y1[2] 0.17 0.00 0.42 0.00 0.00 1.00 1.00 1964 1998
## y1[3] 0.63 0.00 0.79 0.00 0.00 2.00 1.00 1454 1506
## y1[4] 2.05 2.00 1.37 1.48 0.00 5.00 1.00 1802 1841
## y1[5] 4.91 5.00 1.68 1.48 2.00 8.00 1.00 2051 2085
## y1[6] 7.89 8.00 1.37 1.48 5.00 10.00 1.00 1831 NA
## y1[7] 9.33 10.00 0.87 0.00 8.00 10.00 1.00 1490 NA
## y1[8] 9.80 10.00 0.45 0.00 9.00 10.00 1.00 1886 NA
## y1[9] 9.95 10.00 0.23 0.00 9.00 10.00 1.00 1963 NA
## y1[10] 9.99 10.00 0.12 0.00 10.00 10.00 1.00 1935 NA
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=ymed),xy,col='black')
up and down
y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))ex8-6.stan
up and down
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=200> b0;
real<lower=0,upper=1> b1;
real<lower=0,upper=1> b2;
real<lower=0,upper=10> s;
}
model{
y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-6.stan')
fn(mdl,data)## Initial log joint probability = -574.473
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 55 -5.86343 2.19592e-05 0.0029161 1 1 137
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -5.86
## b0 98.77
## b1 0.03
## b2 0.20
## s 0.81
## m0[1] 35.51
## m0[2] 48.72
## m0[3] 34.66
## m0[4] 25.66
## m0[5] 14.61
## m0[6] 37.52
## m0[7] 58.04
## m0[8] 58.46
## m0[9] 25.17
## m0[10] 23.49
## m0[11] 58.10
## m0[12] 60.53
## m0[13] 55.61
## m0[14] 26.04
## m0[15] 53.35
## m0[16] 38.28
## m0[17] 57.80
## m0[18] 23.54
## m0[19] 56.38
## m0[20] 46.36
## y0[1] 36.37
## y0[2] 48.25
## y0[3] 35.77
## y0[4] 26.26
## y0[5] 15.61
## y0[6] 38.06
## y0[7] 59.13
## y0[8] 57.54
## y0[9] 25.62
## y0[10] 22.49
## y0[11] 57.84
## y0[12] 59.27
## y0[13] 53.97
## y0[14] 26.20
## y0[15] 53.73
## y0[16] 38.44
## y0[17] 58.95
## y0[18] 24.36
## y0[19] 57.09
## y0[20] 45.07
## m1[1] 14.61
## m1[2] 54.50
## m1[3] 60.86
## m1[4] 56.89
## m1[5] 50.29
## m1[6] 43.57
## m1[7] 37.45
## m1[8] 32.09
## m1[9] 27.46
## m1[10] 23.49
## y1[1] 14.56
## y1[2] 54.25
## y1[3] 61.37
## y1[4] 56.19
## y1[5] 50.03
## y1[6] 44.03
## y1[7] 36.82
## y1[8] 31.69
## y1[9] 26.92
## y1[10] 24.38
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -9.81 -9.44 1.58 1.36 -12.72 -7.97 1.00 443 565
## b0 98.77 98.73 2.35 2.18 94.93 102.74 1.00 919 550
## b1 0.03 0.03 0.00 0.00 0.03 0.03 1.00 943 692
## b2 0.20 0.20 0.01 0.01 0.19 0.22 1.00 966 662
## s 0.95 0.92 0.18 0.16 0.72 1.31 1.02 256 304
## m0[1] 35.52 35.53 0.28 0.27 35.06 35.99 1.00 1847 1602
## m0[2] 48.74 48.73 0.63 0.61 47.76 49.80 1.00 1245 1317
## m0[3] 34.67 34.67 0.29 0.27 34.20 35.14 1.00 1816 1454
## m0[4] 25.68 25.68 0.35 0.33 25.09 26.26 1.00 1269 1173
## m0[5] 14.62 14.61 0.33 0.32 14.11 15.19 1.00 1060 970
## m0[6] 37.53 37.54 0.28 0.27 37.07 37.99 1.00 1899 1620
## m0[7] 58.04 58.06 0.36 0.34 57.44 58.59 1.00 1264 1081
## m0[8] 58.46 58.48 0.36 0.34 57.87 59.00 1.00 1282 1057
## m0[9] 25.18 25.18 0.35 0.33 24.59 25.77 1.00 1252 1159
## m0[10] 23.51 23.51 0.36 0.34 22.90 24.12 1.00 1210 1028
## m0[11] 58.10 58.10 0.48 0.47 57.35 58.92 1.00 1875 1354
## m0[12] 60.52 60.54 0.35 0.33 59.93 61.06 1.00 1550 1366
## m0[13] 55.60 55.62 0.36 0.34 55.03 56.16 1.00 1222 1080
## m0[14] 26.06 26.06 0.35 0.33 25.48 26.64 1.00 1282 1174
## m0[15] 53.35 53.37 0.35 0.33 52.78 53.90 1.00 1243 1082
## m0[16] 38.29 38.29 0.28 0.27 37.83 38.76 1.00 1906 1593
## m0[17] 57.81 57.81 0.49 0.48 57.04 58.64 1.00 1861 1393
## m0[18] 23.56 23.56 0.36 0.34 22.95 24.17 1.00 1212 1044
## m0[19] 56.37 56.39 0.36 0.34 55.80 56.92 1.00 1228 1092
## m0[20] 46.36 46.37 0.31 0.28 45.84 46.85 1.00 1511 1282
## y0[1] 35.55 35.56 1.02 0.97 33.96 37.24 1.00 1908 1858
## y0[2] 48.72 48.76 1.17 1.11 46.81 50.63 1.00 1801 1776
## y0[3] 34.67 34.70 1.01 0.96 32.97 36.31 1.00 1686 1642
## y0[4] 25.67 25.67 1.02 0.97 24.01 27.35 1.00 1636 1582
## y0[5] 14.60 14.60 1.03 1.02 12.96 16.28 1.00 1876 1704
## y0[6] 37.55 37.55 1.01 0.96 35.90 39.22 1.00 1905 1919
## y0[7] 58.04 58.00 1.05 1.03 56.36 59.72 1.00 1830 1518
## y0[8] 58.48 58.51 1.05 1.00 56.76 60.21 1.00 1853 1555
## y0[9] 25.19 25.21 1.02 0.95 23.49 26.88 1.00 1948 1733
## y0[10] 23.51 23.50 1.03 0.97 21.83 25.16 1.00 1931 1912
## y0[11] 58.12 58.12 1.10 1.02 56.34 59.94 1.00 1703 1808
## y0[12] 60.52 60.53 1.02 0.96 58.82 62.17 1.00 1951 1572
## y0[13] 55.61 55.61 1.06 0.99 53.92 57.37 1.00 1610 1869
## y0[14] 26.04 26.06 1.04 0.98 24.26 27.76 1.00 1974 1729
## y0[15] 53.37 53.37 1.05 0.98 51.70 55.10 1.00 1895 1724
## y0[16] 38.28 38.31 1.04 0.99 36.56 39.90 1.00 2085 1932
## y0[17] 57.78 57.76 1.09 1.06 55.96 59.50 1.00 1950 1739
## y0[18] 23.55 23.53 1.03 1.00 21.88 25.25 1.00 1998 1874
## y0[19] 56.34 56.31 1.04 0.99 54.68 58.07 1.00 1820 1907
## y0[20] 46.35 46.34 1.01 0.97 44.71 48.00 1.00 2046 2057
## m1[1] 14.62 14.61 0.33 0.32 14.11 15.19 1.00 1060 970
## m1[2] 54.51 54.51 0.56 0.54 53.61 55.44 1.00 1542 1378
## m1[3] 60.86 60.87 0.36 0.34 60.26 61.41 1.00 1729 1474
## m1[4] 56.88 56.90 0.36 0.34 56.29 57.43 1.00 1235 1092
## m1[5] 50.29 50.30 0.33 0.32 49.76 50.80 1.00 1319 1271
## m1[6] 43.57 43.57 0.29 0.27 43.08 44.04 1.00 1721 1362
## m1[7] 37.46 37.46 0.28 0.27 36.99 37.92 1.00 1898 1707
## m1[8] 32.10 32.10 0.30 0.28 31.61 32.60 1.00 1603 1354
## m1[9] 27.48 27.48 0.33 0.32 26.92 28.04 1.00 1335 1135
## m1[10] 23.51 23.51 0.36 0.34 22.90 24.12 1.00 1210 1028
## y1[1] 14.61 14.60 1.03 1.00 12.93 16.31 1.00 2009 1645
## y1[2] 54.47 54.50 1.13 1.05 52.57 56.21 1.00 1893 1370
## y1[3] 60.86 60.85 1.06 1.01 59.17 62.61 1.00 2069 1892
## y1[4] 56.89 56.93 1.02 0.99 55.24 58.54 1.00 2006 1823
## y1[5] 50.27 50.26 1.05 1.00 48.52 51.98 1.00 1817 1702
## y1[6] 43.56 43.57 1.01 0.96 41.90 45.20 1.00 1801 1653
## y1[7] 37.47 37.47 1.00 0.94 35.83 39.12 1.00 2123 1770
## y1[8] 32.12 32.12 0.97 0.94 30.55 33.68 1.00 2093 1810
## y1[9] 27.47 27.48 1.02 0.95 25.81 29.16 1.00 1922 1740
## y1[10] 23.48 23.46 1.06 1.04 21.76 25.17 1.00 1491 1530
